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Q ■ Abstract 

u 

We study the influence of diffusion on the scaling properties of the first order struc- 
ture function, Si, of a two-dimensional chaotically advected passive scalar with finite 
lifetime, i.e., with a decaying term in its evolution equation. We obtain an analytical 
^ \ expression for S\ where the dependence on the diffusivity, the decaying coefficient 

o. 

and the stirring due to the chaotic flow is explicitly stated. We show that the presence 
of diffusion introduces a crossover length-scale, the diffusion scale (Lrf), such that 



o 



the scaling behaviour for the structure function is analytical for length-scales shorter 



> 



than and shows a scaling exponent that depends on the decaying term and the 



mixing of the flow for larger scales. Therefore, the scaling exponents for scales 
larger than are not modified with respect to those calculated in the zero diffusion 
limit. Moreover, turns out to be independent of the decaying coeficient, being its 
value the same as for the passive scalar with infinite lifetime. Numerical results sup- 
port our theoretical findings. Our analytical and numerical calculations rest upon the 
Feynmann-Kac representation of the advection-reaction-diffusion partial differential 
equation. 

PACS: 47.52.+j, 05.45.-a, 47.70.Fw, 47.53.+n 
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I. INTRODUCTION 



The subject of advection of (chemically or biologically) reacting substances is of major im- 
portance in fields ranging from combustion to water quality control (see for example [[]]]). One of 
the simplest models of reaction occurring in a flow is the decrease at a constant rate in the concen- 
tration of a substance. Radioactive or photochemical decay belong to this class of reactions. The 
spatial structure of the decaying chemical under chaotic fluid stirring has been described in recent 
works [0,^]] with emphasis on the micro structure formation by the stretching properties of the flow, 
at scales at which molecular diffusion can be neglected. The interpretation of real data, specially in 
laboratory experiments in which the diffusion scale may be explicitely resolved, requires however 
the consideration of diffusive processes which tend to homogeneize the structure at the smallest 
scales. 

In this paper we investigate the role of the molecular diffusion in the scaling properties of the 
first order structure function, defined below (Q) when q = 1, for a passive scalar advected by a 
two-dimensional flow yielding Lagrangian chaos, and with finite lifetime. Thus, the passive scalar 
field (passive in the sense that there is no back influence on the hydrodynamics) 0(x, t) evolves 
according to 

+ v(rj t) . V 0(r, t) = z/V 2 0(r, t) - 60(r, t) + S(v), (1) 

where v is the diffusivity, —b(p represents the decaying and indicates that the passive scalar field 
decays with time at rate b > 0, S(r) is a steady source of the passive scalar, and finally, v(x, t) is 
the velocity of the flow. Here we assume that the velocity field is two-dimensional, incompressible, 
smooth, and nonturbulent. Chaotic advection in two dimensions is obtained generically if a time 
dependence (periodic, for e.g.) is included in v(x, t) |Q]. 

Eq. (jl]) can be considered as an approximation to more complex biological or chemical phe- 
nomena, like the dynamics of plankton populations in ocean currents [|5|^] and the advection of 
pollutants or chemical substances in the atmosphere [0], or as a description of simple procceses, 
like the relaxation of the sea surface temperature [j^] . Moreover, in [g] it is argued that, in a de- 
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termined scale range, the vorticity of a turbulent flow may evolve passively according to Eq. ([]]), 
indicating the high relevance of this equation in the studies of two-dimensional turbulence. 

In the different case of a turbulent velocity field, the spatial propiertes, neglecting intermittency 
corrections, of a scalar field evolving through Eq. ([j]) were first written down by Corrsin [jTI]]. 
The analytical expression for the power spectrum was obtained, among other situations, in the 
Batchelor's regime, that is, between the smallest typical length scale of the velocity field and the 
characteristic length scale of diffusion (La), and from this, the scaling exponents and the diffusion 
length scale, L d , were derived. The calculated form of the power spectrum T(k) is 

r(Jfe) = Nk^ 2b ^- l exp{-2vk 2 /X), (2) 

where N is a constant and A is the absolute value of the most negative of the average stretching 
rates. However, no numerical test of this result, checking the influence of molecular diffusion 
on the crossover, has been reported. Moreover, the criticisms made by Kraichnan JTT| ] to the 
crossover of the power spectrum calculated by Batchelor [jr2p for a passive scalar with infinite 
lifetime, b = 0, as being too sensible to the intermittency, can also be applied here. 

In this paper we study the case of a smooth (nonturbulent) chaotic flow, and instead of the 
power spectrum we calculate the first order structure function. In the absence of diffusion, the 
scaling exponents for the structure functions were calculated neglecting intermittency corrections 
[@], or including them for smooth chaotic flows [|3|JT3|], and in [[14]] for the Kraichnan flow in the 



spatially smooth limit. The structure functions are important quantities than can yield information 
about the typical variation of the concentration field over a small distance Sx. Thus, a structure 
function of order q is defined as 

S q (5x) =< \5<j)\ q >=< |0(x + 5xn) -0(x)| 9 >, (3) 

where <> denotes spatial average, n is a unit vector in the chosen direction, and q is a positive 
real number. In general, for small 5x the structure functions are expected to exhibit a power-law 
dependence S q (Sx) ~ 5x^ q characterized by the set of scaling exponents ( q . Neglecting inter- 
mittency corrections implies that the scaling exponent ( g varies linearly with q. In this work we 



neglect intermittency, and therefore, we can limit ourselves to study only the first order structure 
function, Si(5x). It is important to note that working with the first order structure function will 



minimize the effect of the intermittency [ |T5[ ] when comparing our analytical results with numeri- 
cal (as it is performed here) or real data. In this sense, we improve the calculations performed by 
Corrsin because the intermittency corrections are more important for the power spectrum than for 
the first order structure function. 

Summing up, in this work we study the spatial properties of a scalar field evolving with Eq. 
(|I]) when the velocity field is smooth and chaotic. We calculate the first order structure function 
and check the influence of the diffusion on its scaling exponent and crossover. We also perform 
numerical calculations to check our analytical results. Numerically we show that effectively the 
calculations performed with the first order structure function improve @. Finally, we mention that 
we base all our analytical and numerical calculations in the Feynman-Kac representation of Eq. 
([Q) [ ^6[^7| ]. In other words, we try to obtain information about the spatial properties of an Eulerian 
field 0(x, t) in the long-time limit, through calculations following individual fluid trajectories, i.e., 
using a Lagrangian formulation [JT7p. Our numerical code is explicitly stated in the numerical 
results section. 



II. ANALYTICAL CALCULATIONS OF THE FIRST ORDER STRUCTURE FUNCTION 

The solution of Eq. ([]]), with initial condition </>(x, 0) = 0o( x ) can be written in terms of the 
so-called Feynmann-Kac representation, [ [T^ ] 

0(x, t) = ^ [r(0)]e- w + f Q dse- b ^S[v(s)}^ (4) 

where r(t) is the solution of the Langevin equation 

^ = v(r,t) + V2^T7(t) (5) 

which satisfies the final condition r(£) = x. It is worth to mention that (|]) considers the backwards- 
in-time dynamics and the problem is well-posed by fixing the final conditions instead of the inital 
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ones. rf(t) is a normalized vector-valued white noise term with zero mean, i.e. < rj(t) >= and 
< r)(t)ri(t') >= 15 (t — t'), with I the identity matrix. The average (•)„ is taken over the different 
stochastic trajectories r(i) ending at the stated final point x. Note that 77 is a dummy variable 
in (Q), which disappears after the averaging. In consequence, in expressions containing several 
0's, for example at different space points or times, a different noise variable, each one statistically 
independent from all the others, should be introduced for every appearance of 0. 

It is important to note that in the long-time limit the <p field does not approach a steady dis- 
tribution but one with the same time-dependence of the flow. Thus, for time-periodic flows the 
0(x, t — ► oo) field is also time-periodic. However, its singular characteristics do not change in 
time, that is, it is a statistically steady field [@J3|]. In the following we focus in the situation in 
which the initial concentration is O = 0, so that all the structure arises from the source. Our 
analytical calculations follow closely along the steps of [QJ^J^], thus we submit the reader to these 
references for further details. We proceed by calculating the difference at time t for the values of 
the chemical field at two different points x + 5x/2, and x — <5x/2 separated by a small distance 
<5x. The expression for the difference <50(x, t; 5x) = 0(x + 5x/2, t) — 0(x — 5x/2, t) is 



Here r^s) and r 2 (s), with < s < t, are the solutions of Eq. (Q) with driving noise terms r} 1 and 
rj 2 , respectively, and final values ri(t) = x + 5x/2, and r 2 (t) = x — <5x/2. Two independent noise 
processes rj 1 and rj 2 have been introduced since two concentration values are used in Eq. @. 
We are interested in the asymptotic (t — > oo) scaling of Eq. @ for <5x finite but small. The 
contributions in the integral (||) split into very different behaviors: there is a time t s such that if 
t s < s < t, the backwards trajectories i"i(s) and r 2 (s) exponentially increase its initial distance 
<5x, since they are advected by a chaotic flow. In this regime, the difference 5S will also grow. 
At t s the trajectory separation is of the order of the system size, or of some coherence length of 




(6) 



where 



(7) 
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the advection velocity, L, and thus can not continue to grow. Thus, for s < t s , 5S will stop its 
systematic growing and S'[r 1 (s)] and S'[r 1 (s)] will fluctuate between the range of values taken 
by S in the system, in a manner independent on the initial separation 5x. Since there is no longer 
difference in the statistical properties of Jq" dse~ b ^~ s ^ S'fr^s)] for i — 1,2, the average value of the 
difference of both quantities, which is the contribution to (Q) from the interval (0, t s ) will vanish. 

To further analyze the behavior of the nonvanishing contribution, from (t s , t), we introduce a 
third trajectory r(s) satisfying again the Langevin equation (|5|) with noise ry (independent of rj 1 
and ry 2 ) and endpoint r(t) = x, and consider the time-dependent differences £r-j(s) = r-j(s)— r(s), 
for i = 1, 2 (for < s < t). They satisfy the stochastic equations : 

^i(s) = v(Ti(s),s) - v(r(s), s) + 2v^&(s) • (8) 

The final values are 5i"i(t) = <5x/2 and Sr 2 (t) = — <5x/2, and the new stochastic processes £i(s) = 
(rj^s) — r](s))/ y/2 turn out also to be white noise terms with zero average and correlation matrix 

>= SijIS(s - s'). 

As far as the diferences 8r remain small, the difference in velocity fields in (|S|) may be lin- 
earized so that 

^-6Ti(s) = J(r(s)) ■ 5^(3) + 2v^6(s) • (9) 
as 

J is the Jacobian matrix of the velocity field v. Equation @ can be formally integrated in terms 
of the fundamental matrix M(s, t), s < t, which is the solution of 

^-M(s,t) = J(r(s))-M(M) (10) 
as 

with final condition M(t, t) = I. The result is: 

Ti(s) = r(s) + 6Ti( 8 ) = r(s) ± M(s, t) ■ y + G^s, t) , (11) 
The positive sign applies to i = 1, the negative to i = 2, and 

Gi(s, t) = 2^ J* ds'M(s, s') ■ ti(s') . (12) 



Gj is independent of 5x. If the flow v(x, t) is chaotic, the matrix M will in general produce 
exponential growth of the second term in ( [TT] ) and of the variance of G^. Thus the linearization 
leading to @ will only be justified for \s — t\ small enough. This may be used to define t s as the 
smallest value of s for which linearization is still valid. From the assumed smallness of 5x, and if 
t s < s, the second term in the r.h.s. of ([TT]) will be small, so that one can write 



where again the two signs refer to the two values of i. We are now in conditions to estimate the 
contribution to @ arising from (t s , t): 



Since rf 1 and rj 2 have exactly the same statistical properties, the first two averages in the r.h.s. are 
identical, and cancel out. The next two averages are also identical and add up, so that @, which 
only receives contributions from (t 3 , t), can be written 



The matrix M(s, t) contains the quantitative details on the exponential separation of trajectories. 
In the deterministic dynamics, for which 5r(s) = M(s, t) ■ <5x, the time scale of exponential 
separation is given by the largest Lyapunov exponent A, so that 5x(s) becomes aligned with the 
unit vector along the local unstable direction c[r(s)] in a time of order A -1 , so that 



S[r t (s)]^S[r(s) + G t (s,t)} 

± WS [r(s) + Gi(s, t)} ■ M(s, t)~ + °(^ 2 ) 



(13) 




(14) 



<50(x, t; 5x) 




(15) 
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6r{s) « c[r(s)]e- A(s -* ) c t (x) • 5x (16) 

if \s — t\ > A -1 . c^(x) is the unit vector dual to c(x). Both vectors are time dependent for 
time-dependent flows, but we do not indicate such dependence for notational simplicity. In this 
backward dynamics, the relevant Lyapunov exponent is the most negative one. But in twodimen- 
sional incompressible flow, it is equal in absolute value to the largest positive one in the forward 
dynamics. In ( |T6| ) we take A > and the sign is explicitly written. Analogously, the local ex- 
panding direction in the backwards dynamics is the contracting one in the forward dynamics. 
Expression ( ]T^D will be introduced in (O). It should be said however that the integral in ( jT3| ) 
contains contributions also for \t — s\ < A -1 . The slow convergence of the Lyapunov exponent to 



its asymptotic value Q18D may introduce corrections, specially if b is large. In accordance with our 
aim of neglecting any intermittency corrections, we approximate the action of M by < JT6| ) in all the 
range of integration. 

We now use the mean value theorem to take out of the integral the temporal average of VS • 
c(r(s)), which we call AS. We introduce in < JT5| ) the change of variables u = e A ^ _s ^and write 



t s = t — r (r > 0). On physical grounds r — > oo if 5x and v — > 0, and it is large but finite for 
small 5x and v. Taking the limit t — > oo and regarding that, in the large time limit, the statistical 
properties of the concentration field are steady: 



5(j)(x, oo; <5x) w c(x) t • 5x ( — — f duu *\ (17) 

\\h / r]i 



Since AS is independent of 5x, the scaling with this last quantity is determined by 

<ty(x, oo; <Jx) = Ac(x) t • 5x ((e T{X - b) ) r - l) (18) 



The only stochastic quantity in < J18| ) is r = t — t s , so that the problem has been reduced to the 
calculation of the statistics of this quantity. In terms of the characteristic function of r, defined as 
W(z) = (exp {-zt)) t , Eq. (|18|) reads 

50(x, oo; 5x) = B [W(b - A) - 1] 5x (19) 

where B = Acos(a), with a the angle of 5x with the local expanding direction, and Sx = |5x|. 



Since t s can be estimated from < Jll| ) as |<5r li2 (t s )| ~ the calculation of the statistics of r is a 



classical first passage-time problem, in which one looks for the time it takes a stochastic process 
to reach a given level. There is a vast literature on such problem for stochastic processes of the 
form (JTT]) in which M is a constant matrix or a constant number [jT9|]. 
In the same spirit as before, we use Eq. (jT5|) to approximate ( |TT] ) by 



foils) 



-A(s-t) 



5x 



c[r(s)] ±c(x)+. — + 2^F(s,t) 



(20) 



with 



F(s,t)^J t S ds'c[r(s')]-Us)e 



-A(t-V) 



(21) 



where we have used ( |16| ) inside the integral in (|12|). 

If \t — s\ ^> A -1 , F becomes a constant random number, Foo. From its definition, it is a 
Gaussian random number, of zero average and variance < (-Foo) 2 >= 1/(2A). Thus from ( |20| ) we 
can write the equation defining r: |5rj(i s )| « L, with t s = t — r: 

5x /2// 



C X 



,t 



-5 



2 ' V A 

where g is a random Gaussian number of zero average and unit variance, from which 

1 , 2L 



(22) 



T 



A 



In 



(23) 



We have introduced 7 = cos(a). In the following we denote Ld = 2y which introduces, as 
will be be seen below, the diffusive length scale. The statistics of r may be thus obtained from 
a change of variables from the Gaussian statistics of g. The characteristic function W(z) can be 
now written in terms of (E31) and the Gaussian distribution of g: 

(2L)- 2 / A roo 



W(z) 



271 



f°° -Si 

/ dge 2 \^5x + L d g\ 
J — 00 



(24) 



At this point we mention that expression ( JTPj ) with W given by ( [2"4] ) constitutes the main result 
of this work. 

More explicit expressions can be obtained in the important cases. If Ld <C Sx, then 
W(z) ~ ('jSx I 'L) z / X . In this case, ( |T9| ) scales as Sx 11 , with H = min(l,6/A), in agree- 
ment with previous results On the contrary, for small scales, that is, if 5x L d , then 



W(z) = 7r 1 / 2 (4i//L 2 A) z / A r(^tA), and H — 1, thus 50 is smooth below a diffusive scale given 
by L d = Therefore, we realise of the existence of a crossover length-scale, the diffu- 



sive scale, given by L d which separates a diffusion-controled smooth behaviour of 5(p from an 
advection-controled with scaling exponent given by min(l, b/X). This last is non-smooth when 
b< A. 

We can obtain an approximate expression simpler than (|23|) by realizing that the fluctuations 
in r are much smaller than the average value if Sx and v are small, so that < e' ZT >~ e~ z<T> , 
and then substitute the average value < t > by an estimation of it, r, obtained for example by the 
condition that the second moment of the trajectory separation reaches the system size: 



,2 At 



„2Ar 



(25) 



from which 



t = — In 
A 



2L 



which is an approximation to the mean first passage time < r > if 5x and u/X are small. Thus 



(26) 



and 



W{z) 



x, oo; 5x) ~ C cos(a) 



7 2 5x 2 + L\ 
~2L 



(7 2 / 2 + 11) 



(27) 



b-\ 
2A 



b-X 

2~ 



I 



(28) 



We have introduced I = Sx/L and l d = L d /L, and C is a constant equals to AL/2~ . 

The first order structure function, Si, can be obtained by averaging (^8j) over the different 
spatial points x. Here we realise that the point spatial dependence of ( |28| ) is only through the angle 
a between <5x and the local expanding direction. Assuming that these directions are isotropically 
distributed, Si is calculated as: 



Si(l) =< \5<P(xl, oo; 5x)\ > 



r2 Tf 




/ da 


cos(a) 


Jo 





(cos(a) 2 / 2 + 1\ 
10 



b-\ 

2\ 



2~ 



(29) 



The integral can be performed and we obtain: 

Si(l) = —I ( (f 2 + if) 2A F (±, |; 

+ 2^), (30) 

being F(a, b, c; x) the confluent hypergeometric function [|Z0"|]. The important thing to be noted 
in ( ffl3| ) is that shows the same scaling behaviour as S(p, with a crossover length-scale given 
by l d = 2^f/L. That is, St(l) ~ Z for / << Z d and Si(i) ~ ^ when / >> l d and b < \. 
Remarkably, the diffusive scale, l d , is independent of the decaying coeficient, b, and therefore its 
value is the same to the one obtained for the passive scalar with infinite lifetime, i.e., with 6 = 0. 
Next section is devoted to check numerically some of the above analytical results. 

III. NUMERICAL RESULTS 

In this section we will check the two regimes found for the scaling of the structure function, 
and also the expression Eq. (j3Q|). 

Numerically, we proceed by integrating backwards in time Eq. (|5]) with initial conditions on a 
one-dimensional transect of the scalar field. Many different trajectories are obtained in this way 
for different realizations of the noise. Then we integrate Eq. (Q) forward in time for each one 
of these trajectories (0 O = 0), and finally, we average over the different trajectories. With this 
procedure we obtain the field </>(x, t), at a fixed time t large enough, on a one-dimensional transect 
without having to calculate the whole two-dimensional field. Obviously, this allows us to reach a 
high resolution for the structure functions which are then only calculated on a ld transect. For the 
flow, v = (v x , v y ), we take a simple time-periodic velocity field defined, in the unit square, L = 1, 

and with periodic boundary conditions, by 

2U T 

v x (x,y,t) = — Y®(~2 ~ t m ° d T ) cos ( 27r ^) 

v y (x, y, t) = —^-&(t mod T — —J cos(27rx) (31) 

where Q(x) is the Heavy side step function, and T is the periodicity of the flow. In our simulations 
U = 1.0, which produces a flow with a single connected chaotic region in the advection dynamics. 
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The value of the numerically obtained Lyapunov exponent is A rs 2.35/T. The source term used 
is S(x,y) = 0.2 sin(27rx) sin(27n/). We perform our calculations until a final typical time of 
tf ~ 15T where we realise that a final statistical stationary concentration field is reached. 

In Fig. (1) we show (in logarithmic scale) the structure functions against the length-scale 
in units of the l d . We plot the curves obtained for different values of the decaying coefficient 
b = 0.5,0.9, 1.2, and for a fixed value of the diffusivity v = 5 x 10~ 7 , and T = 1. For any 
of the different plots we observe three different regimes. First, for small scales and up to the 
diffusion length scale l d (\n(l/l d ) = 0) the slope of the curve is 1, showing that the structure 
function is smooth in this length-scale interval. Then, for scales larger than l d and up to typical 
scales comparable to the system size L — 1, appears the other scaling regime. A straight line with 
slope | is observed here. In the figure, we also plot the straight lines showing this behaviour, and 
with the number above indicating their slope. 
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_7 1 i i i i i i i i i i 

-2 2 4 6 8 

ln(l/l d ) 

FIG. 1. This figure plots ln(Si(Z)) against \n(l/l d ). Always T = 1 and v = 5 x 1CT 7 . Here b = 1.2 
for the dotted line, b = 0.9 for the dashed line, and b = 0.5 for the long-dashed line. The straight line 
above each curve shows the scaling behaviour, and the number over the line indicates its slope. The 
two regimes are observed: a) smooth (slope 1 for \n(l/l d ) < 0, b) scaling behaviour with slope b/X for 
ln(l) < ln(l/l d ) « \n(L/L d ) = hx(l/l d ). 

Next, we proceed to check numerically the validity of Eq.(^). First, we realize that, on di- 
mensional grounds, the adimensional ratio Si = S±(l) / Si(l = l d ) should be a function just of the 
adimensional parameters I, l d , and b/X. Expression ( PTj| ) fully satisfies this requirement. This can 
be used to compare numerically-obtained values of Si as a function of I (or 1/ l d as in the plot) 
for a variety of values of A, b, and v with the theoretical prediction. This last function would be 
a single curve as long as the combinations l d = 2^2v/X/L and b/X do not change. Figure (2) 
shows (for two different values of (u/X,b/X)) that the numerical data indeed scale as expected 
from the theory, and that the data collapse into the analytic curves is very good, confirming the 
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accuracy of expression (|30|). 

Finally, in Fig. (3) we perform an analogous plot (in logarithmic scale) to Fig. (2a) but 
for the power spectrum (also normalised to the diffusion length scale) analytically calculated by 
Corrsin (0), and the ones obtained numerically for the same values of b, A and v. Just a visual 
inspection confirms our expectation that, due to the effect of intermittency, it is more reliable to 
compare with numerical or real data the approximations to the first order structure function than 
the approximations to the power spectrum. Also, it is worth to mention that to calculate the power 
spectrum much more data are needed to obtain a nice statistics than in the case of the first order 
structure function. 




-4 -2 2 4 6 -4 -2 2 4 6 
ln(l/l d ) ln(l/l d ) 
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FIG. 2. Plot of the adimensional structure function, Si, as a function of I /Id = Sx/Ld- With solid line 
we plot the theoretical result (Eq.(|30|)), and with different symbols the different values obtained numeri- 
cally, a) is for the fixed values b/X = 0.382 and u/X = 4.2 x 1CT 6 . b) the same for b/X = 0.639 and 
v j X = 2. 1 x 10~ 5 . The symbols in the legends indicate the values of b, X and v used to obtain the numerical 
data. 
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FIG. 3. Comparison of the power spectrum (solid line) obtained by Corrsin (|J) and numerical data as 
in Fig. (2a). 

IV. SUMMARY 



To summarize, in this paper we have investigated the role of diffusion in the scaling properties 
of the first order structure function of a passive scalar chaotically advected. Based on the Feynman- 
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Kac representation of the advection-reaction-diffusion equation, we have obtained an analytical 
expression for the the shape of the first order structure function for all scales smaller than the 
typical system size, and observed the existence of a crossover between a smooth regime and a 
scaling regime with slope depending on the rate of the decaying coefficient and the maximum 
Lyapunov exponent of the flow. The length-scale for the crossover is the typical dissipative length 
scale, given by L d = 2 which is independent of the decaying coeficient and, therefore, is 
the same as for the passive scalar with infinite lifetime. Moreover, our results confirm that the 
scaling exponent of S\ are not modified by the molecular diffusion. We have shown that our 
approximation for the first order structure function is more reliable than the one for the power 
spectrum obtained under the same disregard of intermittency corrections. This could be relevant 
when explaining or modeling numerical or experimental data. We have also provided numerical 
support to our analytical findings. These are in excellent agreement. Our numerical algorithm 
uses the Feynman-Kac representation, and allows for a very fine scale resolution, as we just need 
to calculate the scalar field in one-dimensional sections of the whole surface. We believe that this 
numerical algorithm can be very useful for other kind of advection-reaction-diffusion problems. 
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